
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/copy.h>
#include <thrust/fill.h>
#include <thrust/sequence.h>
#include <thrust/transform_reduce.h>
#include <thrust/functional.h>


__global__ void insertKernel(int *a)
{
	int x = blockDim.x * blockIdx.x + threadIdx.x;
	a[x] = x;
	
}


__global__ void initKernel(int *a, int *b, double penal, double *Di, int column, int laa, int lbb)
{
	int x = blockDim.x * blockIdx.x + threadIdx.x;
	int i,j;
	i = x % lbb;
	j = x / lbb;
	if ((i < laa) && (j < lbb))
	{
		if (a[i] == b[j])
		{
			Di[((i + 1)*column) + j + 1] = 1;
		}
		else
		{
			Di[((i + 1)*column) + j + 1] = -penal;
		}
	}
}
__global__ void verticalscanKernel(double *acc_cst, int column, double Lh, double Lv, double penal, int roww, double *temp_maxx, int ii, int jj, int *bestpredii,int *bestpredijj)
{
	int x = blockDim.x * blockIdx.x + threadIdx.x;
	
	if ((x >= roww) && (x<ii - 1))
		{
			if ((acc_cst[(x*column) + jj]) - (1 + (penal)*(ii - x)) > *temp_maxx)
			{
				*temp_maxx = (acc_cst[(x*column) + jj]) - (1 + (penal)*(ii - x)); //the second term is the penalty term for the vertical transition
				*bestpredii = x;
				*bestpredijj = jj;
			}
		}
}
__global__ void horizontalscalscanKernel(double *acc_cst, int column, double Lh, double Lv, double penal, int roww,int coluu, double *temp_maxx, int ii, int jj, int *bestpredii, int *bestpredijj)
{
	int x = blockDim.x * blockIdx.x + threadIdx.x;

	if ((x >= coluu) && (x<jj - 1))
	{
		if ((acc_cst[(ii*column) + x]) - (1 + (penal)*(jj - x)) > *temp_maxx)
		{
			*temp_maxx = (acc_cst[(ii*column) + x]) - (1 + (penal)*(jj - x)); //the second term is the penalty term for the vertical transition
			*bestpredii = ii;
			*bestpredijj = x;
		}
	}
}

__global__ void gridKernel(double *Di,double *acc_cst, double *row_pre, double *col_pre, int column, double Lh, double Lv, double penal, int roww)
{
	int x = blockDim.x * blockIdx.x + threadIdx.x;
	int i, j, bestPredi, bestPredj,rowstart, row, col, colstart;
	double temp_max;
	i = x % column;
	j = x / column;
	if ((i >= 1) && (j >= 1) && (i < roww) && (j < column))
	{
		temp_max = Di[(i*column) + j];
		bestPredi = 0;
		bestPredj = 0;
		//Diagonal transition
		if ((acc_cst[((i - 1)*column) + j - 1] + Di[(i*column) + j]) > temp_max)
		{
			temp_max = acc_cst[((i - 1)*column) + j - 1] + Di[(i*column) + j];
			bestPredi = i - 1;
			bestPredj = j - 1;
		}
		//Vertical scan: nodes (1,j),(2,j),...,(i-1,j)
		if ((i - (int)Lv) >= 1)
			rowstart = i - Lv;
		else rowstart = 1;
		for (row = rowstart; row < i - 1; row++)
		{
			if ((acc_cst[(row*column) + j]) - (1 + (penal)*(i - row)) > temp_max)
			{
				temp_max = (acc_cst[(row*column) + j]) - (1 + (penal)*(i - row)); //the second term is the penalty term for the vertical transition
				bestPredi = row;
				bestPredj = j;
			}
		}
		//Horizontal scan: nodes (i,1),(i,2),...,(i,j-1)
		if ((j - (int)Lh) >= 1)
			colstart = j - Lh;
		else colstart = 1;
		for (col = colstart; col < j - 1; col++)
		{
			if ((acc_cst[(i*column) + col]) - (1 + (penal)*(j - col)) > temp_max)
			{
				temp_max = (acc_cst[(i*column) + col]) - (1 + (penal)*(j - col)); //the second term is the penalty term for the vertical transition
				bestPredi = i;
				bestPredj = col;
			}
		}
		//Finished (i,j).There only remains to store the winner
		if (temp_max > 0)
		{
			acc_cst[(i*column) + j] = temp_max;
			row_pre[(i*column) + j] = bestPredi;
			col_pre[(i*column) + j] = bestPredj;
		}
	}
}

int main()
{
	int counter1, patterncounter, i, LA, M, LB, N, size, size2, x, j, bestPredi, bestPredj,rowstart,colstart,row,col;
	int data_to_read;
	FILE *data1; FILE *pattern;
	int *temp_memory; int *temp_pattern; int *device_vector1; int *device_pattern;
	double penalty, Lhor, Lver,maxv,temp_max;
	double *H; double *D; double *temp_acc_cost; double *acc_cost; double *temp_row_pred; double *row_pred;
	double *temp_col_pred; double *col_pred;
	//number initiallizations
	counter1 = 0; patterncounter = 0;
	i = 0;
	maxv = 0.0;

	clock_t tic = clock();

	/*files opening and reading size*/
	if ((data1 = fopen("data1.txt", "r")) == NULL)
	{
		fprintf(stderr, "error in opening file:%s", "data1.txt");
		exit(-1);
	}
	while (!feof(data1))
	{
		fscanf(data1, "%d", &data_to_read);
		counter1++;
	}
	fclose(data1);
	if ((pattern = fopen("pattern.txt", "r")) == NULL)
	{
		fprintf(stderr, "error in opening file:%s", "pattern.txt");
		exit(-1);
	}
	while (!feof(pattern))
	{
		fscanf(pattern, "%d", &data_to_read);
		patterncounter++;
	}
	fclose(pattern);
	// temp memory vector has storage for counter1 integers
	temp_memory = (int*)malloc(counter1*sizeof(int));
	temp_pattern = (int*)malloc(patterncounter*sizeof(int));
	std::cout << "temp_memory has size " << counter1 << std::endl;
	std::cout << "temp_pattern has size " << patterncounter << std::endl;
	//filling temp_memory with systemcall id
	if ((data1 = fopen("data1.txt", "r")) == NULL)
	{
		fprintf(stderr, "error in opening file:%s", "data1.txt");
		exit(-1);
	}
	i = 0;
	while (!feof(data1))
	{
		fscanf(data1, "%d", &temp_memory[i]);
		//std::cout <<temp_memory[i] << std::endl;
		i++;
	}
	fclose(data1);
	if ((pattern = fopen("pattern.txt", "r")) == NULL)
	{
		fprintf(stderr, "error in opening file:%s", "pattern.txt");
		exit(-1);
	}
	i = 0;
	while (!feof(pattern))
	{
		fscanf(pattern, "%d", &temp_pattern[i]);
		i++;
	}
	fclose(pattern);
	//device_vector1 for pattern and data_set1
	cudaMalloc((void**)&device_vector1, counter1*sizeof(int));
	cudaMemcpy(device_vector1, temp_memory, counter1*sizeof(int), cudaMemcpyHostToDevice);
	cudaMalloc((void**)&device_pattern, (int)patterncounter*sizeof(int));
	cudaMemcpy(device_pattern, temp_pattern, patterncounter*sizeof(int), cudaMemcpyHostToDevice);
		
	//smith waterman start
	Lver = 5.0;
	Lhor = 5.0;
	penalty = 1.0 / 3.0;
	LA = counter1;
	M = LA + 1;
	LB = patterncounter;
	N = LB + 1;
	size = N*M;
	size2 = LA*LB;
	H = (double*)malloc(size*sizeof(double));
	cudaMalloc((void**)&D, size*sizeof(double));
	for (i = 0; i < size; i++)
	{
		H[i] = 0;
	}
	cudaMemcpy(D, H, N*M*sizeof(double), cudaMemcpyHostToDevice);
	// print contents of D 
	printf("%f \n", penalty);
	//system("pause");
	int threadsPerBlock = 1024;
	int blocksPerGrid = (size2 + threadsPerBlock - 1) / threadsPerBlock;
	initKernel <<<blocksPerGrid, threadsPerBlock >>>(device_vector1, device_pattern, penalty, D, N, LA, LB);
	cudaMemcpy(H, D, size*sizeof(double), cudaMemcpyDeviceToHost);
		
	temp_acc_cost = (double*)malloc(size*sizeof(double));
	cudaMalloc((void**)&acc_cost, size*sizeof(double));
	for (i = 0; i < size; i++)
	{
		temp_acc_cost[i] = 0;
	}
	cudaMemcpy(acc_cost, temp_acc_cost, size*sizeof(double), cudaMemcpyHostToDevice);

	temp_row_pred = (double*)malloc(size*sizeof(double));
	cudaMalloc((void**)&row_pred, size*sizeof(double));
	for (i = 0; i < size; i++)
	{
		temp_row_pred[i] = 0;
	}
	cudaMemcpy(row_pred, temp_row_pred, size*sizeof(double), cudaMemcpyHostToDevice);

	temp_col_pred = (double*)malloc(size*sizeof(double));
	cudaMalloc((void**)&col_pred, size*sizeof(double));
	for (i = 0; i < size; i++)
	{
		temp_col_pred[i] = 0;
	}
	cudaMemcpy(col_pred, temp_col_pred, size*sizeof(double), cudaMemcpyHostToDevice);
	
	//start grid processing
	for (x = 0; x < size;x++)
	{
		i = x % N;
		j = x / N;
		if ((i >= 1) && (j >= 1) && (i < M) && (j < N))
		{
			temp_max = H[(i*N) + j];
			bestPredi = 0;
			bestPredj = 0;
			//Diagonal transition
			if ((temp_acc_cost[((i - 1)*N) + j - 1] + H[(i*N) + j]) > temp_max)
			{
				temp_max = temp_acc_cost[((i - 1)*N) + j - 1] + H[(i*N) + j];
				bestPredi = i - 1;
				bestPredj = j - 1;
			}
			
			//Vertical scan: nodes (1,j),(2,j),...,(i-1,j)
			if ((i - (int)Lver) >= 1)
				rowstart = i - Lver;
			else rowstart = 1;
			cudaMemcpy(acc_cost, temp_acc_cost, size*sizeof(double), cudaMemcpyHostToDevice);
			cudaMemcpy(row_pred, temp_row_pred, size*sizeof(double), cudaMemcpyHostToDevice);
			cudaMemcpy(col_pred, temp_col_pred, size*sizeof(double), cudaMemcpyHostToDevice);
			blocksPerGrid = (size + threadsPerBlock - 1) / threadsPerBlock;
			verticalscanKernel << <blocksPerGrid, threadsPerBlock >> >(acc_cost, N, Lhor, Lver, penalty, rowstart, &temp_max, i, j, &bestPredi, &bestPredj);
			cudaMemcpy(temp_acc_cost, acc_cost, size*sizeof(double), cudaMemcpyDeviceToHost);
			cudaMemcpy(temp_row_pred, row_pred, size*sizeof(double), cudaMemcpyDeviceToHost);
			cudaMemcpy(temp_col_pred, col_pred, size*sizeof(double), cudaMemcpyDeviceToHost);
						
			//Horizontal scan: nodes (i,1),(i,2),...,(i,j-1)
			if ((j - (int)Lhor) >= 1)
				colstart = j - Lhor;
			else colstart = 1;
			cudaMemcpy(acc_cost, temp_acc_cost, size*sizeof(double), cudaMemcpyHostToDevice);
			cudaMemcpy(row_pred, temp_row_pred, size*sizeof(double), cudaMemcpyHostToDevice);
			cudaMemcpy(col_pred, temp_col_pred, size*sizeof(double), cudaMemcpyHostToDevice);
			blocksPerGrid = (size + threadsPerBlock - 1) / threadsPerBlock;
			horizontalscalscanKernel << <blocksPerGrid, threadsPerBlock >> >(acc_cost, N, Lhor, Lver, penalty, rowstart, colstart, &temp_max, i, j, &bestPredi, &bestPredj);
			cudaMemcpy(temp_acc_cost, acc_cost, size*sizeof(double), cudaMemcpyDeviceToHost);
			cudaMemcpy(temp_row_pred, row_pred, size*sizeof(double), cudaMemcpyDeviceToHost);
			cudaMemcpy(temp_col_pred, col_pred, size*sizeof(double), cudaMemcpyDeviceToHost);
						
			//Finished (i,j).There only remains to store the winner
			if (temp_max > 0)
			{
				temp_acc_cost[(i*N) + j] = temp_max;
				temp_row_pred[(i*N) + j] = bestPredi;
				temp_col_pred[(i*N) + j] = bestPredj;
			}
		}
	}
	//eo grid processing
	maxv = temp_acc_cost[0];
	for (i = 1; i < size; i++)
	{
		if (temp_acc_cost[i] > maxv)
			maxv = temp_acc_cost[i];
	}
	
	printf("the similarity is %f \n", maxv);
	clock_t toc = clock();

	printf("Elapsed: %f seconds\n", (double)(toc - tic) / CLOCKS_PER_SEC);


	system("pause");


   	// empty the vector
	cudaFree(device_vector1);
	cudaFree(device_pattern);
	free(temp_memory);
	free(temp_pattern);
	cudaFree(device_vector1);
	cudaFree(device_pattern);
	cudaFree(acc_cost);
	cudaFree(col_pred);
	cudaFree(row_pred);
	free(temp_acc_cost);
	free(temp_col_pred);
	free(temp_row_pred);
		return 0;
}